@read_strat
;-------------------------------------------------
; basic parameters for base state density profile
; must agree with the values taken in the eulag
; version used for the corresponding simulation 
;-------------------------------------------------
; add:
; 
; they are naturally configured to the eulag 
; values by defining them with the initial
; values of strat.dat
; --> still need to reset the value of "g"
;_________________________________________________

rh00=rho_strat(0)
t00=temp_strat(0)
th00=t00
;g=1.
if g eq !Null then begin
 print, 'please set the value of the g aceleration ...'
endif else begin
 rg=0.4
 cp=2.5*rg
 cap=rg/cp
 capi=1./cap
 ss=g/(cp*th00)
 press00=rg*rh00*t00
 gam=5./3.

;-------------------------------------------------
; Radial discretization -> must agree with 
; /Downloads/code/rk4PACHO.f of the eulag version
; used for the corresponding simulation.
;-------------------------------------------------

 zmin=z_strat(0)
 zmax=z_strat(-1)
 dz=(zmax-zmin)/l
 z=indgen(l)*dz

;-------------------------------------------------
;     ADIABATIC PROFILES !!! 
;
; --> Base state density, pressure, temperature
;     & and theta profiles.
;
; --> Local adiabatic speed of sound profile
;
; --> notice that in this case the theta profile
;     reduces to the case of constant entropy:
;
;     theta_ad = t_ad*((t00*rh00)/(t_ad*rho_ad))^cap
;              = t00    
;
;-------------------------------------------------

 rho_ad = rh00*(1-ss*(z))^(capi-1)
 p_ad = press00*(1-ss*z)^(capi)
 t_ad = t00*(1-ss*z)
 ;theta_ad = t_ad*((t00*rh00)/(t_ad*rho_ad))^cap
 theta_ad = t00
 cs_ad=sqrt(gam*p_ad/rho_ad)
  
;________________________________________________________________
; SPATIAL AVERAGES -> gives variables of the type f=f(z) only,
;                     denoted as <f>
;________________________________________________________________ 
; --> FOR THE ENERGY FLUXES:
;________________________________________________________________
;  The convective and kinetic fluxes are given by:
;
;                      F_kin = (0.5)*(rho_strat)*<u^2*w>
;                      F_conv = (cp)*(rho_strat)*<w*DT>
;
;  Where the horizontally averaged temperature fluctuations are 
;  given by:
;                      DT = (T_ad/theta_ad)*D(theta)
; 
;  Where theta is the total potential temperature and D(theta) are
;  the corresponding fluctuations identified with the symbol Th'
;  in the EULAG literature.
;________________________________________________________________
; --> FOR MIXING LENGTH:
;________________________________________________________________
;
;  The reference comparison to be made with MLT is given by:
;
;          DT/T ~ <w^2>/cs^2 ~ [F_conv/(rho_strat*cs^3)]^2/3 
;
;  The relative temperature fluctuations are given by:
;
;                      DT/T = <Th'>/(<Th'>+theta_strat)
;
;  the speed of sound is given by the usual definition:
;
;                      cs=sqrt(dP/d(rho))
;
;  which for this case is given by the polytropic atmosphere
;  and so we get:
;  
;                      cs=sqrt((1+1/n)*(P_strat)/(rho_strat))
;________________________________________________________________  

; --> arrays of the type array[z,t], with z the radial coordinate
;     and t being the time

 the_aver =total(the,1)/m
 w2_aver  =total(w2,1)/m
 thew_aver=total(thew,1)/m
 p_aver   =total(p,1)/m
 ;u2w_aver =total(u2w,1)/m      ;for the kin flux (still not implemented)
                                 ;in eulag

 print, 'horizontally averaged variables ( w[z,t] ) computed:'
 print, 'the_aver -> <Th>'
 print, 'w2_aver -> <w2>'
 print, 'thew_aver -> <w*Th>'
 print, 'p_aver -> <p>'
;-------------------------------------------------
; Computes horizontally averaged fluxes
;------------------------------------------------- 
 
 z_num=n_elements(the_aver(*,0))
 t_num=n_elements(the_aver(0,*))
 
 fconv=make_array(z_num,t_num)
 fkin=make_array(z_num,t_num)
 
 for ii=0,t_num-1 do begin
  for jj=0,z_num-1 do begin
   fconv(jj,ii)=cp*rho_strat(jj)*(t_ad(jj)/theta_ad)*thew_aver(jj,ii)
   ;fkin(jj,ii)=0.5*rho_strat(jj)*u2w_aver(jj,ii)
  endfor
 endfor 
 print, 'fconv -> Fc = cp*rho_e*(To/Tho)*<w*Th>'
 print, '________________________________________________________________'

endelse  
